On the unsteady behavior of turbulence models 
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Periodically forced turbulence is used as a test-case to evaluate the predictions of two-equation 
and multiple-scale turbulence models in unsteady flows. The limitations of the two-equation model 
are shown to originate in the basic assumption of spectral equilibrium. A multiple-scale model 
based on a picture of stepwise energy cascade overcomes some of these limitations, but the absence 
of nonlocal interactions proves to lead to poor predictions of the time variation of the dissipation 
rate. A new multiple-scale model that includes nonlocal interactions is proposed and shown to 
reproduce the main features of the frequency response correctly. 



A basic premise of one point closures such as the k — e 
model is the hypothesis of "spectral equilibrium," which 
justifies two distinct roles of the dissipation rate e: on 
the one hand, it appears in the energy balance, defined 
as a correlation of velocity gradients, hence a small-scale 
quantity; on the other hand, it is used phenomenolog- 
ically to describe large-scale transport properties. The 
most basic formulation of the latter is Kolmogorov's hy- 
pothesis e oc fc 3 / 2 /L where k is the turbulent kinetic 
energy and where I is a length scale characteristic of 
the largest scales of motion; equivalent formulations in- 
clude v t oc fc 2 /e, the formula for turbulent viscosity, or 
t oc k/e, the formula for the turbulent time-scale. The 
Kolmogorov theory, or more general assumptions of self- 
similarity of all scales of motion, justify all of these pro- 
portionalities although the constants of proportion- 
ality need not coincide in all self-similar flows [2] . 

But turbulence models are not needed to describe self- 
similar flows, which merely serve as calibration cases; 
models are needed to describe departure from self- 
similarity, when spectral equilibrium becomes a strong 
constraint on turbulence evolution. In a study of a flow 
in which turbulence evolves from a steady state to a self- 
similar time-dependent state Q, the implications of the 
departure from spectral equilibrium were investigated: 
this departure was connected to transient failure of the 
Tennekes-Lumley balance Q and to the consequent rel- 
evance of small scale dynamics for the large scales. The 
breakdown of spectral equilibrium has also been observed 
in engineering flows including turbulent diffusers and 
wakes. [1] 

The limitations of the spectral equilibrium hypothe- 
sis are very well known in the modeling literature and 
have led to proposals for multiple- scale models Q that 
more realistically address the complexity of the non- 
linear interactions in turbulent flows. This letter re- 
ports on some investigations of multiple-scale models ap- 
plied to an especially simple and attractive test case for 
transient turbulence: periodically forced turbulence . 
This problem arises when isotropic incompressible tur- 
bulence, maintained in a steady state by a large-scale 
isotropic forcing with total amplitude p, is subjected to 
a small time-dependent periodic perturbation with am- 
plitude p: p{t) = p + p cos(wi), such that the ratio 



p/p <C 1 and such that the forcing length scale does 
not depend on time. The phase-averaged kinetic en- 
ergy k can then be decomposed into a mean k and a 
periodic part fc(oj) cos(wt + 0fc(u;)), with <fik the phase- 
shift between the forcing and the kinetic energy. Sim- 
ilarly, the viscous dissipation rate can be written as 
e = I + e(uj) cos(W + c (w)); k and e are related to the 
time-independent forcing length scale L by L oc k 3 ^ 2 /e. 
The periodic parts of k and e are sinusoidal, like the forc- 
ing, because p/p <C 1. The functions k(cu), e(w), 0fc(w), 
and 4>e{oj), which characterize the linearized response of 
steady-state turbulence to periodic perturbation of the 
forcing, can be called the linear response junctions. 

We will use periodically forced turbulence as a test case 
to evaluate the ability of multiple-scale models to predict 
the dynamics of time-dependent turbulence. It will first 
be shown that the unsteady predictions of multiple-scale 
models are significantly better than the predictions of a 
two-equation model. However, an elementary multiple- 
scale model based on the heuristic picture of stepwise 
energy cascade is found to have limitations in predicting 
the unsteady dissipation rate. A multiple-scale model 
that includes the possibility of nonlocal interactions is 
proposed; it is shown that this model can capture some 
fine features of the unsteady energy dissipation. 

The linear response functions were determined in re- 
cent work [lfj using the EDQNM closure theory, which 
was shown to compare very well to available low Reynolds 
number experimental, and DNS data. Comparison with 
high Reynolds number data for linear response func- 
tions would be desirable, but such data is not yet avail- 
able. Briefly summarizing the major conclusions, the 
two-equation model is satisfactory both in the static limit 
u) — > 0, in which the phase shifts 4>k, (fre vanish, and in the 
frozen limit u> — > oc, in which k ~ uj^ 1 and 4>k ~ w/2. 
However, this agreement is trivial; only the results at in- 
termediate frequencies provide a real test of the model. 
The two-equation model reproduces the function k(u>) 
reasonably well, but the transition from the static limit 
(f>k ~ to the frozen limit <fik ~ 7r/2 occurs over a fre- 
quency range that is much too wide. The amplitude e(w) 



is not satisfactory: a range in which e 



at high 



Reynolds numbers is absent. We note that since observa- 
tions of the modulated dissipation rate are very difficult, 
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for the time derivative of the spectral energy flux F(k, t), 



FIG. 1: Schematic view of a discretized energy cascade. The 
solid arrows indicate the energy fluxes between neighboring 
wavenumber shells. The dashed lines indicate the nonlocal 
fluxes between spectrally remote shells. Note that a similar 
picture can be found in a paper by Lumley (Tlj . 



and relevant high Reynolds number data is not yet avail- 
able, the EDQNM results for this quantity remain the- 
oretical predictions; they are nevertheless supported by 
arguments based on the well-established role of dis- 
tant interactions in turbulence. Finally, the two-equation 
model also makes the incorrect prediction that 4>k — 4>e- 
regardless of u. 

We investigate whether these predictions can be im- 
proved using multiple-scale modeling following the ideas 
of Schiestel [6J. Multiple-scale models can be considered 
numerical methods for spectral closures, with significant 
modifications designed to permit reasonable accuracy at 
a very low order of discretization. Thus, whereas a nu- 
merical implementation of a spectral closure would solve 
for the energy spectrum at perhaps hundreds of discrete 
wavenumbers «,-, in Schiestel's formulation, the energy 
spectrum is divided into a relatively small number of 
wavenumber shells Ki-i < k, < Kf, for each shell, equa- 
tions are written for two scalar descriptors: the fluctu- 
ation energy contained in the shell and the net energy 
flux into it. To enhance the accuracy possible with a 
relatively small number of shells, Schiestel allowed the 
partition wavenumbers m to be functions of time. A 
schematic picture of the resulting discretized energy cas- 
cade is given in Figure [JJ following Lumley [Tll j . 

The starting point for the analytical formulation is the 
Lin equation governing the energy spectrum E{n,t), 



d_ 
dt 



— + vk 2 ) E(n,t) = P(K,t) 



dF(K,t) 
dn 



(1) 



In this equation v is the viscosity, P(k, t) the forcing 
term, and F (k, t) is the energy flux across wavenumber 
K. Schiestel applied the Kovaznay model [l2f 



F(«,i) = CK 5/2 E{ K ,t) 3/2 , 



(2) 



with C a parameter which determines the Kolmogorov 
constant. This model represents a stepwise cascade of 
energy in spectral space from small to large k and repro- 
duces a k -5 / 3 inertial range. Figure [T] would correspond 
to a stepwise cascade if the dashed lines were absent. 
To obtain a multiple-scale model, write the equation 



F(K,t) 



3F(n,t) 
2 E(k, t) 



E(n,t). 



(3) 



From the viewpoint of Schiestel's analysis, we have as- 
sumed that the partition wavenumbers are constant in 
time: this assumption seems appropriate for this prob- 
lem, in which the forcing wavenumber is fixed. The pro- 
duction scales are assumed to be confined to low n, and 
the dissipation to high k (a more general description 
including finite-Reynolds number effects or broad-band 
forcing will not be attempted here). For n in the inertial 
range, ((T|) becomes 



E{k, t) 



dF(K,t) 
8k 



We combine this with ^ and © to obtain 

V ' ; 2E(n,t) dn 



(4) 



(5) 



A discrete model is obtained by splitting the spectral do- 
main into n shells as illustrated in Figure [TJ The energy 
in shell i is e% « £7(/tj)A/ti. The spectral flux /, ss F(ni) 
and the time derivative is dF(ni,t)/dk w (/j— /i_i)/A/tj, 
so that we obtain 



fi 



n (fi fi—lj 

2 e. 



1< i < n. 



(6) 



Integrating ^ over each shell gives the partial energy 
balance equations 



-(fi-fi-l), l<i<n. 



(7) 



in which all fi and ei are functions only of time. The os- 
cillating production term p, assumed to act at the small 
wavenumbers, is identified with /o, the flux entering shell 
I. Furthermore, the high Reynolds numbers case is con- 
sidered in which we assume that the viscous dissipation 
takes place at the last wavenumber shell: e = f n . This as- 
sumption makes it unnecessary to introduce partial dissi- 
pation rates for each shell and corresponding equations of 
motion. A special feature of the Kovaznay model is that 
the partition wavenumbers do not appear in the model. 

By choosing n = 1, one obtains a two-equation model; 
the equation for fi becomes the dissipation rate equa- 
tion e = (3/2)(e/fc)(P - e). Note that the two model 
constants, generally called C e i and C e 2 in the literature, 
are equal, which allows the study of statistically station- 
ary isotropic turbulence. Consistency with homogeneous 
shear flow, or with any problem in which the forcing 
length scale increases as a power law or exponential in 
time [III, requires C e i < C e 2- 

Choosing n > 2 should improve the predictions by in- 
troducing the possibility of spectral imbalance, a nec- 
essary requirement if the same model is to be applied 
to both forced and decaying turbulence Q; imbalance is 
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possible because the partial fluxes fi with 1 < i < n — 1 
and the dissipation e = f n are independent. 

The unsteady behavior of this model is now assessed 
as follows. Beginning with a steady state with shell en- 
ergies e-i such that p = /j = e, the periodic perturbation 
pcos(wi) is added to the forcing. Describing the periodic 
response in terms of complex amplitudes k(to) and e(w), 
so that k(uj) = \k(uj)\ and tan0 fe (w) = 3fc(cj)/3?A;(w) 
with the obvious analogs for e, the equations for the pe- 
riodic part of the partial kinetic energies and local fluxes 
become 

3 f. „ 

iuie.i = -(/i- fi-i), iuifi = -~(fi - fi-i). (8) 

2 ej 

The resulting linear system is easily solved analytically 
for ii(ui) and fi(uj) in terms of the periodic forcing per- 
turbation p and the parameters and fi . 

Figure [5] compares the results for the amplitude k(u) 
and phase shift 4>k{^->) for models with n — 1 to n = 7 
wavenumber shells. Also shown are results obtained in 
EDQNM computations [Io| at Reynolds number R\ = 
1000, with R x « 15i? L 1/2 and i? L = (2k/3) 1 / 2 L/v. It 
has been shown [l(| that at low u, k(u)), should tend 
to a plateau and that at large uj, k(uj) follows a power- 
law proportional to lu" 1 . Confirming the conclusion 
of Bos et aZ.[T3], k(u>) is in reasonable agreement with 
EDQNM even for the two-equation model n = 1, al- 
though this agreement improves significantly as the num- 
ber of wavenumber partitions n increases. 

The error in the two-equation model prediction of k(u>) 
is a too gradual transition from the plateau to the lu^ 1 re- 
gion. This defect appears more prominently in the phase 
</>fc(u): although all models give the correct static and 
frozen limits, the two-equation model transitions much 
too gradually, and only the models with n > 2 are in close 
agreement with EDQNM. The relatively rapid transition 
in the energy phase shift therefore appears as a typical 
multiple-scale effect: evidently, in this problem, the small 
scales are not simply 'slaved' to the large scales through a 
constant dissipation rate as is assumed in a two-equation 
model; instead, they are dynamically independent and 
have a strong effect on what is apparently a purely large- 
scale property. 

It has also been demonstrated [Io[ that the response 
function e(uj) follows at high uj a power-law proportional 
to u~ 3 up to the Kolmogorov frequency uj v ~ e/v. For 
uj > ujjj, e{ui) becomes proportional to These results 
are shown for comparison in the graph on the left side 
of Figure [3] (the phase shift 4> e proves difficult to com- 
pute with any confidence because of the extremely small 
amplitudes involved, therefore comparisons are omitted) . 

The agreement of e(w) given by the multiple-scale 
model Eq. © with the EDQNM results is very good 
down to values e/p = 10~ 3 for n > 3 (note that k(uj) and 
e(uj) are proportional to p, and p is chosen unity without 
loss of generality). However for smaller values of e/p, the 
discrete model starts to diverge from the EDQNM re- 
sults, especially for large n. It is very easily shown that 
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FIG. 2: Linear response functions k(ui) (left) and </>fc(aj) 
(right) for the multiple-scale model based on the Kovaznay 
spectral closure ©. The 7 curves in each graph correspond 
to models with n = 1, 2, .., 7 spectral shells and n increases in 
the direction of the arrow. The theoretical amplitude predic- 
tion uj^ 1 is indicated by a straight line. Also shown are the 
results of EDQNM simulations (symbols). The frequency uj 
is normalized by the large scale frequency e/k. 



for this model, the leading order contributions at high oj 
are proportional to ui~ n : indeed, recursive solution of the 
equations for fi in ([8]) show that fi ~ ui~ l p. Thus, only 
if n = 3 can we obtain £(oj) ~ w~ 3 , but this is the result 
of a coincidence, which disappears if the number of par- 
tition wavenumbers is increased. This difficulty reflects 
a limitation of the multiple-scale model: Eq. implies 
a linear first-order partial differential equation for F in 
which disturbances in F propagate along characteristics; 
this property is probably significantly compromised by a 
finite dimensional approximate model. 

It has been shown[l(| that nonlocal interactions are 
responsible for the range. Nonlocal interactions are 
represented in Figure[T]by dashed lines. Such interactions 
do not occur in the model Eqs. ([6|) and (O, because 
quantities in shell i depend only on its nearest neighbor, 
shell i — 1. The absence of nonlocal interactions will be 
even more significant in problems in which the role of 
nonlocality is greater, as in the Batchelor regime of the 
passive scalar [lj] and in MHD [15j . 

To address this problem, a new multiple-scale model 
will now be derived including the effect of nonlocal inter- 
actions. We start from a simple spectral model contain- 
ing nonlocal interactions, due to Ellison 

F(K,t) = c(J K 2 E{n,t)df?j nE{n,t). (9) 

Applying the same procedure as to the Kovaznay model, 
the partial energy balance Eq. Q is unchanged, but Eq. 
(IH) is replaced by 



e,: 



fi J2 P =i K l(fp fp-i) 

~ J2 P =i K p e p 



(10) 



This expression contains the wavenumbers /tj because 
nonlocal interactions, which depend on the spacing be- 
tween the wave-number partitions, have been retained. 
The specification of the Ki becomes part of the model: 
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FIG. 3: Frequency response function e(w) for the model based 
on the Kovaznay closure © (left) and for the model based on 
the Ellison closure (right). The 7 curves in each graph 
correspond to n = 1,2, ..,7 and n increases in the direction 
of the arrow. The theoretical prediction u)~ 3 is indicated by 
the straight lines. Also shown are the results of EDQNM 
simulations (symbols). 

we will use a logarithmic discretization with 
in which r is a model parameter which determines the 
logarithmic gridsize. Using this discretization, the ra- 
tio K n /m — r n ~~ l , so that a large range of scales can 
be considered by increasing r. Note that if a linear dis- 
cretization is used, K n j n\ — (tiA/c/IAk) = n so that the 
number of partitions for high Reynolds numbers becomes 
prohibitively large. Using the logarithmic discretization, 
the model for periodic forcing becomes 

* _ fi ( i t x fi Y, P =i r2{p 1 KU - fp-i) , , 

e-i I T, P =i r(p > e P 

with the same partial energy equations as in Eq. (|8]l. 
Again, when n = 1 the model reduces to a two-equation 
model. For n > 1 the model differs from the previous 
model through the interaction term which couples wave- 
number shell i with all wavenumber shells p = 1, ■ • • , i. 
The response functions depend on both n and r, which 



will be chosen through a compromise between computa- 
tional cost and precision. 

The results for e(ui) obtained from Eq. (fTTj) with r = 
3 and 1 < n < 7 wavenumber partitions are shown in 
the graph on the right side of Figure [3] In the same 
figure, EDQNM results (lf| at R\ = 1000 are shown. 
The agreement with EDQNM is very good down to values 
e/p = 10~ 3 for 3 < n < 7, and for n = 7, agreement is 
good down to e/p — 10~ 6 . We conclude that the model 
including nonlocal interactions Eq. (|lip makes better 
predictions of e than the model Eq. ([§])■ in which nonlocal 
interactions are absent. The predictions of Eq. (jlll) (not 
shown) for fc(w) and </>/c(w) very nearly coincide with the 
results obtained using Eq. ([5]). 

To conclude, we have found that in the problem of 
periodically forced turbulence, a two-equation model 
only gives satisfactory predictions for k(oj) and <fik(w) at 
asymptotically high and low frequencies. The predic- 
tions at intermediate frequencies are improved by using 
a multiple-scale model based on the heuristic picture of 
stepwise energy cascade, the Kovaznay model. In par- 
ticular, this model correctly predicts the rapid jump of 
the phase shift </>/c(w) between the static and frozen lim- 
its. The multiple-scale model based on the Ellison clo- 
sure includes nonlocal effects, and leads to better agree- 
ment with EDQNM, including the high Reynolds num- 
ber w~ 3 scaling range [Toj. Both models give practically 
indistinguishable results for the amplitude and phase of 
the oscillating kinetic energy, which is not strongly influ- 
enced by nonlocal interactions. Our approach suggests 
how one might construct reduced order models for phe- 
nomena dominated by significant nonlocal interactions, 
like the Batchelor range of the passive scalar and some 
cases of MHD. 

We would like to acknowledge the interesting and 
thoughtful comments of the referees, which led to sig- 
nificant modifications of the paper. 
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